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Abstract 


A technique to obtain the sensitivity of the static aeroelastic response of a three di- 
mensional wing model is designed and implemented. The formulation is quite general 
and accepts any aerodynamic and structural analysis capability. A program to combine 
the discipline level, or local, sensitivities into global sensitivity derivatives is developed. 
A variety of representations of the wing pressure field are developed and tested to de- 
termine the most accurate and efficient scheme for representing the field outside of the 
aerodvnamic code. Chebyshev polynomials are used to globally fit the pressure field. This 
approach had some difficulties in representing local variations in the field, so a variety 
of local interpolation polynomial pressure representations are also implemented. These 
panel based representations use a constant pressure value, a bilinearly interpolated value, 
or a biquadratically interpolated value. The interpolation polynomial approaches do an 
excellent job of reducing the numerical problems of the global approach for comparable 
computational effort. Regardless of the pressure representation used, sensitivity and re- 
sponse results with excellent accuracy have been produced for large integrated quantities 
such as wing tip deflection and trim angle of attack. The sensitivities of such things as 
individual generalized displacements have been found with fair accuracy. In general, ac- 
curacy is found to be proportional to the relative size of the derivatives to the quantity 
itself. 
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1.0 Introduction 


During the design phase of an engineering system, numerous analyses are conducted to 
predict changes in the characteristics of the system due to changes in design variables. Usu- 
ally, this process entails perturbing each variable in turn, recalculating the characteristics, 
and evaluating the sensitivities by a finite-difference calculation. These repeated analyses 
can drive the cost of design very high. An approach that has found increased interest 
recently in engineering design is analytical calculation of the sensitivity derivatives 1 . Typ- 
ically, the analytical approach requires less computational resources than a finite-difference 
approach and is less subject to numerical errors (round-off or truncation). The analytical 
approach is best developed in parallel with the baseline analysis capability since it uses a 
significant portion of the numerical information generated during that baseline analysis. 
In the design of modern aircraft, airframe flexibility is a concern from strength, control, 
and performance standpoints. To properly account for the aerodynamic and structural im- 
plications of flexibility, reliable aeroelastic sensitivity analysis is needed. Therefore, both 
structural and aerodynamic sensitivity analysis capabilities are necessary. 


1.1 Sensitivity Analysis 


Structural sensitivity analysis methodology has been available for over two decades 
for both sizing (thickness, cross-section properties) and shape (configuration) variables . 
However, aerodynamic sensitivity analysis has been nonexistent until relatively recently. 
Some limited aerodynamic sensitivity analysis capability was developed for aircraft in 
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subcriticai compressible flow by Hawk and Bristow 3 , but it only handled perturbations in 
the direction of the thickness of the wing (thickness, camber, or twist distribution). Yates 4 
proposed a new approach that considers general geometry variations including planform 
for subsonic, sonic, and supersonic unsteady, nonplanar lifting-surface theory. This work 

is still under progress. 

Aeroelastic sensitivity analysis methodology has also been available for more than 
two decades for structural sizing variables (Haftka and Yates 5 ). This is because changes 
in sizing variables exclusively affect the structural stiffness and mass distribution of the 
airframe and not its basic geometry. Therefore, structural sensitivity analysis capability 
is sufficient. However, the lack of development in aerodynamic shape sensitivity analysis 
explains why there are very few results in aeroelastic shape sensitivity analysis. In a 
notable exception. Haftka et al 6 designed a sailplane wing under aeroelastic constraints 
and analyzed the design model with vortex lattice and finite element methods. A finite- 
difference (for GSE derivatives) aeroelastic sensitivity analysis capability is made possible 
by (1) devising a reduced order model to describe the wing static aeroelastic response and 
(2) using exact perturbation analysis to approximate changes in the \ortieity \ector with 
changes in the geometry. Follow up work used semi-analytical derivatives. 

Barthelemy and Bergen 7 demonstrated the feasibility of calculating analytically the 
sensitivity of wing static aeroelastic characteristics to changes in wing shape. Of interest 
also was the fact that the curvature of the aeroelastic characteristics was small enough 
that analytical sensitivity derivatives could be used to approximate them without costly 
reanalyses for large perturbations of the design variables. 

The dynamic aeroelastic phenomena is also of interest to designers and it would be 
advantageous to the aircraft designers to have a tool that can be used to predict the changes 
in flutter speed with the changes in basic shape parameters. 

As is the case for static aeroelastic response, sensitivity calculations have only been 
available for structural sizing parameters. For example, Rudisill and Bhatia 3 developed 
expressions for the analytical derivatives of the eigenvalues, reduced frequency and flutter 
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speed with respect to structural parameters for use in minimizing the total mass. However, 
this method is limited because the structural parameters are sizing variables such as cross- 
sectional areas, plate thickness and diameters of spars. 

Pedersen and Seyranian 9 examined the change in flutter load as a function of change 
in stiffness, mass, boundary conditions or load distribution. They showed how sensitivity 
analysis can be performed without any new eigenvalue analysis. The solutton to the mam 
and adjoint problem prov.de all the necessary information for evaluating sensitivit.es. Their 
paper mainly focuses on column and beam critical load distributions. 

In a recent study, Kapania. Bergen and Barthelemy 19 obtained the sensitivity of 
a wing flutter response to changes in its geometry. Specifically, the objective was to 
determine the derivatives of flutter speed and frequency with respect to wing area, aspect 
ratio, taper ratio, and sweep angle. The study used Giles-'- equivalent plate model to 
represent the wing structure. The aerodynamic loads were obtained using Yates' 3 mod.fied 
strip analysis to analyze flutter characteristics for finite span swept and unswept wings. 

It is noted that Yates modified strip theory was used quite recently by Landsberger and 
Dugundji 14 , with a mod, Scat, on for camber effects g,ven by Spielberg 15 , to study the flutter 
and divergence of a composite plate. 

Unger, Hutchison, Rais-Rohan., Haftka. and Grossman 15 recently demonstrated a 
variable-complexity approach to the mult.disciplinary design of a transport wing. Their 
approach optimizes based on a relatively simple and inexpensive model. Then, this ap- 
proximate optimal wing is reoptimized using a more complex model. 

An excellent examination of various issues involved in aeroelast.c analysis was re- 
cently published by Borland 17 . The study compares and contrasts “integrated” aeroelast.c 
systems versus “interfaced” systems. Borland refers to packages des.gned to produce aeroe- 
lastic results from a single program as “integrated”. Packages that have separate modules 
for each discipline that each write as their output the input for another package are termed 
■interfaced". As is discussed in section 2.1 of this dissertation, the current work's scheme 
is a variation on the “interfaced” approach. 


1.0 Introduction 


1.2 Current Work 


The current work develops techniques to determine the sensitivity of the various static 
aeroelastic responses to the variations in various shape parameters, namely: (i) wing area, 
(ii) sweep, (iii) aspect ratio and (iv) the taper ratio. The aeroelastic responses were the 
generalized aeroelastic displacements and the trim angle of attack. The sensitivities were 
obtained by differentiating the constitutive equations. It was shown that the resulting sen- 
sitivity equations can be reformulated into a variation of the Sobieski’s Global Sensitivity 
Equations 19 (G.S.E.) approach. Both schemes gave the various global sensitivities (i.e. 
the sensitivity including all interdisciplinary interactions) in terms of local sensitivities 
(i.e. the sensitivities obtained at the discipline level). A key feature that distinguished 
this study from the study by Barthelemy and Bergen 7 is the use of a more realistic aerody- 
namic model, FAST 20 , that uses a lifting surface theory as opposed to a lifting line theory 
employed in the earlier study. The formulation was designed to be quite general so that 
it was applicable with any aerodynamic code which, for a given geometry and structural 
deformations, provides aerodynamic pressures on the wing surface. To facilitate the calcu- 
lation of the shape sensitivities of various quantities (required in aeroelastic analyses), the 
pressure distribution was represented as a double series of global Chebyshev polynomials. 
The displacements of the wing were obtained using an iterative scheme. To validate this 
more general formulation, sensitivity of the static aeroelastic response of an example wing 
was obtained. The results were compared with those obtained by using a purely finite 
difference approach. A good agreement was obtained. 

During the course of the work [Ref. 18], it was found that the generalized pressure 
coefficients, due to the global nature of the interpolation polynomials, may be sensitive 
to small changes in independent variables. As a result, the determination of the local 
derivatives of some of the generalized aerodynamic coefficients was found to be difficult 
when forward or central differences were used. A higher order finite difference scheme using 
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a large step size was emplyed so that the effect of local wiggles can be reduced. Obviously, 
this was an expensive option. Other more robust techniques are therefore needed to express 
the aerodynamic pressure. A representation of the aerodynamic loads that is piece-wise 
polynomial should be used. Some of the advantages of piece-wise representation over global 

representation are discussed by Burden and Faires 21 and by de Boor 22 . 

In addition, the piece- wise polynomial pressure approach will have a major application 
in representing the aerodynic pressure of the High Speed Civil Transport Wing. This wing's 
flow regime will be transonic, thereby having pressure discontinuities from shock waves. 
A piece-wise function is often better represented by piece-wise polynomials as opposed to 
globally defined polynomials. 
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2.0 Aeroelastic Response 


A scheme was developed to interface arbitrary aerodynamic and structural codes 
in order to calculate the aeroelastic response of a wing. Anticipating the availability 
of nonlinear aerodynamic models in the future, the formulation does not assume a linear 
dependence between the lift generated and the generalized coordinates and the initial angle 
of attack. This formulation needs an iterative process to calculate the angle of attack to 

which the aircraft is trimmed to produce the required lift. 

This analysis is performed for the wing shown in figure 1. This forward swept wing 
has its coordinate origin at the root quarter chord. The y-axis is perpendicular to the root 
chord, in the span-wise direction. The x-axis is in the chord-wise direction, with positive 
values being aft of the quarter-chord. 


2.1 Overview 


The combination of realistic aerodynamics and structural models in a modular man- 
ner with shape sensitivity code requires a systematic approach. A scheme of calling the 
aerodynamic and structural codes to produce a converged static wing loading and shape 
was developed. In addition a set of “neutral format" data files were defined. Here, a 
“neutral format” data file is a file that is defined to contain certain data at certain spots, 
regardless of the package that originally generated the data. This scheme makes the re- 
placement of analysis packages practical and relatively simple. These files include the base 
geometry and initial deflection values, the intermediate pressure loading and structural 
deflections, and the final converged wing loading and deflection. 
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The aeroelastic problem is broken into subproblems (or blocks) by discipline. The 
aerodynamic and structural blocks are called iteratively to produce a converged static wing 
loading and shape. Shape sensitivity values for this converged wing are then obtained 
(discussed in chapter 3). Figure 1 is a flow chart illustrating this scheme. 

Each block operates completely independently of the other. It reads from several 
“neutral format” input files, performs its calculations and generates one or more neutral 
format output files. Thus, any aerodynamic and structural analysis capability may be 
used. Only new input and output “adapter” programs need be written to add a new 
analysis package to the system. These two adapter programs must convert the neutral 
format data files to and from the new package’s native format. 

2.11 Aerodynamics: 

The aerodynamic block is responsible for generating the loads on the wing. It reads as 
input the wing geometry parameters and the current wing deflections. It is able to output 
the pressure on the wing at arbitrary points. Figure 2 illustrates the aerodynamic section. 
Currently, the aerodynamic analysis is being performed by program FAST. This lifting 
surface code was developed at NASA, Langley. It is based on a kernal-function based 
theory developed by Yates, and was implemented by Desmarais and Bennett 20 . Originally 
developed for a CDC Cyber computer, this program, with considerable effort, has been 
ported to the VAX/ VMS and IBM- CMS operating systems. Adapter programs to convert 
to and from native FAST data files to the defined neutral format were also developed. 

2.12 Structures: 

The structural block is responsible for calculating the deflection of the wing. It is 
given the wing geometry and wing loading. It calculates the deflected shape of the wing. 
Figure 3 illustrates the structural section. Currently, Giles' ELAPS 12 code is being used 
to perform the structural analysis. This Ritz method program was developed at NASA. 
Langley. It has been adapted for use on both the VAX and IBM systems. Adapter 
programs have been developed to convert wing pressures to ELAPS generalized forces and 
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Figure 2. Overall Flow Chart 
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Figure 3. Aerodynamic Section 
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to convert its deflection outputs to neutral form. 


2.2 Governing Equations 


The governing equations of motion for the aeroelastic anal\ sis and the lift can be 


written as 


(A-]{C} = {0} 


(2.1) 


and 


~T~ = J j^p( x *y) d & 


where 

[A'] is the stiffness matrix. 

{C} is the vector of unknown generalized displacements. 
{Q} is the vector of generalized forces. 
n is the load factor, 

W is the aircraft weight, 

p(x,y ) is the wing pressure field, and 

ft is the wing surface area. 


The vector of generalized forces can be obtained as: 


Qi = 



p( x -y)~/i( x 1 y)dzdy 


( 2 . 3 ) 


where 7i(x,y) is the ith nondimensional displacement function used in the displacement 
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model: 


np 

’-''.Ji y) = y. -fi(x,y)Ci (2.4) 

1=1 

with 

= (—)"•( — r (2.0) 

'Emax ymax 

where m varies from 0 to mw and n varies from 2 to nw 4 - 2 with n varying more rapidly 
them m. The terms mw and nw are configurable options. For this work, mw = 5 and 
nw = 6 . 

These 7 ’s satisfy the geometric boundary conditions for a cantilever plate. The C s 
in (2.4) are the generalized displacements. 

To facilitate both the integration and subsequent sensitivity calculations, a coordinate 
transformation was used to simplify the integration limits. This was accomplished using 
the following transformation: 


J (»7.s) = (~-6) 

;=i 

4 

y{n<£) = (2.7) 

j=i 

where 

Nj(t hO are the shape functions, and the 

Xj and y ; are the coordinates of the four corner points of the wing. 

The shape functions are given as 

N,(r } ,0 = ( l + «i)(l+W)/4 ( 2 . 8 ) 

where q, and are the coordinates of the node i in the q — £ system. Xote that this 
transformation will change the domain of the wing to a square (— 1 < 17 < 1 ; — 1 < £ < 1 ). 
The transformed wing is shown in figure 5. 
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As a first step to obtain the generalized forces, the pressure distribution on the wing 
was represented as a series of either global interpolation functions or interpolation poly- 
nomials that have local support over a wing panel. This can be represented in a generic 
form as: 


w 

p(n,Z) = (2.9) 

;= i 

where a- 7 can be considered as the generalized pressure coefficients and j3 J (r},£) are some 
known interpolation functions of r/ and if. 

A large number of interpolating polynomials are available in the literature 21,22 . First, 
a global tensor product of Chebyshev polynomials was used. Then, a variety of local 
polynomials were studied. Details of these approaches and their results are in chapters 
4-7. 

Generically. the integral for a generalized force Q,. is 

Qi = j J PitfiO 7i( r 7-f) |«A( ; ?-s)l drjd£ (2.10) 

where |Ji(r/,£)| th e Jacobian of the coordinate transformation specified by Eqs. 2.6 and 

2.7. The generalized force Q x can be written as 


M 

Q. = Z A ->°’ 

;=i 


In matrix form 

{Q} = [A}{a} 

Similarly, the lift equation can be written as: 


nW 

2 


M 

}= 1 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 
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2.3 Aeroelastic Response 


The aeroelastic response was obtained in an iterative fashion. In that, the pressure 
distribution on the wing is first obtained by assuming the wing to be rigid and having an 
angle of attack of 1° (throughout the span). The pressure distribution thus obtained is 
used to obtain the vector of generalized forces (Eq. 2.9) which in turn is used to obtain the 
vector of generalized displacements (Eq. 2.1). The elastic displacements are superimposed 
on the rigid wing and a new pressure distribution on the wing is obtained. This pressure 
distribution is then used to obtain the generalized displacements. The total lift on the wing 
is calculated, and a new trim angle of attack is obtained by dividing the total required lift 
by the current calculated lift and multiplying by the current trim angle of attack. This 
process is repeated until a converged value of the trim angle of attack is achieved for the 
wing. No relaxation is necessary to achieve convergence for the cases studied. 

An alternative technique has also been developed to compute the trim angle of attack 
and wing generalized deflections. This technique involves significantly more setup effort, 
but gives the converged displacement and angle of attack in one step. It should be noted 
that only the linearity of the two current analysis programs makes this one step solution 

possible. 

For the case of linear aerodynamics, one is not required to use an iterative scheme to 
obtain the aeroelastic solution. For such a case the vector of aerodynamic pressures. p„ at 
some discrete points cam be written as 

( P} = [iHC} + o {|^} (2 ’ 14) 


where = 


dp, 

dc v 


The generalized aerodynamic coefficients can be written as 


{a} = [R] M 


(2.15) 
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where R is an interpolation matrix that converts discrete p’s to our generalized a's. For 
our case, {p} is taken at Chebychev points. Eqn. 2.13 is substituted into eqn. 2.14. Then, 
the equilibrium equations for the structure are updated from eqns. 2.1 and 2.8. 


[A-]{C} = (.4] (i?l [-4] (C) +<*[.■»] [*]{f£} 


The trim equation is adapted from eqn. 2.12. 


= (L} r [fl] [i] fC} + a {in*] 


Note that for our case: [R\ [-4] = 


d{ Q } 

5{C} 


The governing equations, for the case of linear aerodynamics, are: 


(2.16) 


(2.17) 


rrn-H] M 

(■•!] ' 


f C l 

f 0 1 







^ < 


3^ 




i <* j 

nlV 
V 2 


(2. IS) 


The left hand side matrix is identical to the sensivity matrix to be derived in the next 
chapter. If linear aerodynamics are to be used and a sensitivity analysis is desired, this 
result can result in a slight savings of computational effort over the iterative approach. 

In order to get a feel for the performance of the aeroelastic system, the spanload was 
plotted at a couple of constant chord positions for both rigid and flexible wings. Figure 
6 shows the spanload along the nodes closest to the leading edge of the wing. Figure i 
shows the spanload along a spanwise line approximately along the wing quarter chord. The 
differences between the flexible and rigid wings are small, indicating a rather stiff wing. 
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Figure 6. Pressure near Wing Leading Edge 
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3.0 Senstitivity Equations 


A scheme for determining the global sensitivities of the aeroelastxc system to changes 
in shape parameters is developed. Tins scheme can be shown to be a variation of Sobieski s 
Global Sensitivity Equations. A variety of local sensitivity data is combined to produce 
global sensitivity results. Here, the term local sensitivity refers to the sensitivity of an item 
within a particular discipline, such as the sensitivity of the wing stiffness matrix to a change 
in wing sweep. Global sensitivities are dependent on the interaction of the disciplines. A 
global sensitivity example is the variation of wing deflected shape with respect to a change 
in wing taper ratio. 

The goal of this analysis is to produce values for the global sensitivities dC/dr t and 
da/dr t . Here, r, represents any of the wing shape parameters: area, aspect ratio, sweep, 
or taper ratio. Equations 2.1, 2.9, and 2.12 can be used to perform the shape sensitivity 
analysis of static aeroelastic response. In the following development d()/d{) indicates a 
local, single discipline, term and d()/d() indicates a global or total derivative. Taking 
derivatives of the equilibrium and the trim equation, with respect to the shape variable r, 
(namely sweep, aspect ratio, wing area, taper ratio), we obtain 



(3.1) 


(3.2) 


Note that the derivative of lift is zero, because we require the total lift acting on the 
wing to remain constant. Additionally, L is purely a function of geometry, so its partial 
derivative is the same as its total derivative. 
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The vector {dQ/dri} can be obtained as 


dQi_ 

dri 


M 


= £([ 


dA tj 

dr i 


1 4 U \ 

a 3 + Aij — ) 


da 3 

dri 


(3.3) 


where 

da 3 da 3 ^ da 3 dC n da J da 

dri dri + <9C n dr t + da dri 

n — 1 


(3.4) 


and where 

5a J /dr/ is the local sensitivity of the aerodynamic generalized pressure coefficients 
and can be obtained while performing the aerodynamic analysis; 

da 3 jdC n is the derivative of the generalized pressure coefficient with respect to a 
generalized displacement C„, and 

da /dr / is the derivative of the trim angle of attack with respect to r/. 

Note that the matrix [A] is also purely a function of geometry: thus its partial and 
total derivative are the same. 

In matrix form, the global sensitivity of generalized forces becomes 


. dQ . . dA , , . da u d{a} ,,dC da da 

- W w + [4 > W + [ - 41! aTcT H ] + ^'' 4|( 


(3.5) 


It is noted that the major computational expense is the determination of d{a)/d{C). 

The sensitivity of generalized displacements, therefore, becomes (combining Eqs. 3.1 
and 3.5) 


w - w®)] <£> - 't> + [a]{ £, ] + £ lA]{i £ } - 0 !c) 


(3.6) 


In this equation, all the terms on the right hand side are known except for da /dri. 
This can be obtained by considering the sensitivity of the lift equation. Eq. 3.2. This 
equation can be written as: 


, dL i t f , , r i T < da , r nTr l r dC . uuc . . ~p . . 

W <“ } + {L)T{-} + {L )Tl^J-\{-) + — {£} {-} 


da 

dri 


da 

da 


(3.T) 
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The required sensitivity derivatives can be then obtained by simultaneously solving 
the sets of eqns 3.6 and 3.7. With a little rearranging they become: 


[A'] - [A][J$f^] 

<£} r lf{^] 


H - { ^) T {a)-{L) T {^) I 


dri 
da 
dr t ) 


Generating this left hand “sensitivity” matrix by finite differences is extremely corn- 
putationally expensive. The term 9{a)/S{C} is particularly expensive. It reqmres calling 
the aerodynamic code one time for each element in the {C} vector. The sensitiv.ty matrix 
is valid for a particular base geometry, regardless of which r, is of interest. It is only 
generated once and saved for future uses. The left hand sensit.vrtv matrix may be used 
with a different right hand side to determine the converged wing aeroelastic response. Tins 

equation is given by Eqn 2.18. 




12 . 18 ] 


If a sensitivity calculation is to be performed as well, then this approach to deter- 
mining the aeroelastic response is more efficient. If only the converged wing configu ration 
is needed, the iterative approach is slightly better. The iterative approach is also more 
general; it can be directly applied to non-linear aerodynamic codes, whereas Eqn. 2. IS » 
only applicable for linear aerodynamics. 


3.1 Companion Between the Present Formulation and Sobzeskt’s Global Sen 
sitivities Equations (GSE) 


In Ref. 18, Sobieski presented two different formulations (GSE1 and GSE2) to obtain 
the global sensitivities of a multi-disciplinary system in terms of the sensitivtt.es of the 
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subsystems, called local sensitivities. It is of interest to compare the present formulation 
(Eq. 3.8) with Sobieski's Global Sensitivity Equations. 

In terms of the Sobieski's first formulation, called GSE1, the governing equations for 

the wing can be written as: 

Trim: 


fa = L T {{r ( },a) {a} - 


nW 

2 


0 


( 3 . 10 ) 


Aerodynamics: 


{/. (a, W. {C})> = {«} 


(3.11) 


Structures: 


{fc} = [A’(r,)]{C} - [A(r,)]{a} = {0} 

Then the GSE’s are: 


df„ 

da 

dfa 

d{a) 

AL- 1 

5{C} 


f da > 
dr t 


dja 

dr t 

fa) 

da 

a{a} 

d{h) 

d{C] 

< 

itel 

dr t 

> = - < 

d±h± 

drt 

iilcl 

da 

d\fc] 

d{a) 

d{fc } 
d{C) . 


ii£l 

< drt J 


djfcl 

dr t 


■ 0 

{L) T 

o - 


( A 

dr t 



a{h] 

da 

a\h) 

d{a} 

9{fa) 

d{C) 

< 

d {a) 

drt 

> = - < 

<£L1 > 

dr t 

. 0 

-[A] 

[K\ - 


d{ C) 

' drt ' 




Note that, since {/„} is not directly accessable, the above form of the Sobieski s Global 
Sensitivity Equations could not be used. 

In terms of Sobieski’s GSE2 formulation, the governing equations, for the system at 


hand are: 
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Trim: 


(3.15) 


a = a ({a}, r t ) 

Aerodynamics: 

fa) = {a(a, {C}, r,)} (3.16) 

Structures: 

(C) = [AT{r,)]-‘ (A(r,)] fa) (3.1T) 

The global sensitivity equations are: 


■ / 

da 

da 1 


f da ' 


f 3°* > 

d{a) 

d{C) 


dr t 


dr t 

9{a) 

da 

I 

d{a] 

d{C} 

< 

ihl 

dr ( 

► = < 

d { a i . 

dri 

cu 

ft) 

1 

1 

21£ 1 
9{a) 

I 


d{ C) 

dri > 


d{C) 

dr ( 7 


Using Eqs. 3.15 - 3.17, the above set of equations becomes 


■ / 

da 0 1 


( da y 


' 0 

d{a) U 


dr t 


_2iiL 

da 

T d( a ) 

1 d{c) 

< 

ihl 

dr t 

) - < 

d .i a i . 

dri > 

. 0 

i 

i 


d\C) 

< drt * 




Since da/d{a } cannot be found easily, this formulation could not be used. The first 
equation from 3.14 and the last two equations from 3.19 can be combinted to form a system 
as follows: 


a{°) 

8-1 

8 (g) 

3r, 



M -[Ml 

-[A']' 1 [A] [/] 

-{L} t 0 


da 


[ 0 ] 

0 



(3.20) 


For the current study, the global derivatives of {a} aren’t necessary, so eliminating 
the first line of 3.20 and rearranging it will give us equation 3.9. Thus, our formulation is 
a mix of GSE1 and GSE2. 
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3.2 Solving the System 


The terms [A'],[A],{a},{C}, and {1} are all quantities known from the converged base- 
line configuration. A finite difference technique is used to calculate the terms , d{a}/da 
d[K]/dr h and d{a}/dr,. The terms d[A\/dr t and d{L}/dn are computed analytically 
The terms d{C)/dr l and da/dr t are the unknowns found by solving this system of equa- 
tions. 


3.3 Derivatives of [A] and {L} 


These derivatives are calculated by using the chain rule and the derivatives of each 
of their elements. Each element's derivative will be reduced to derivatives of one of our 
coordinates (u, tj, etc.). Thus, derivatives of these coordinates will be found first, then the 

derivatives of the elements will be found. 

Coordinate Derivatives 

u,u - Local panel coordinates. The panel based pressure interpolations discussed 
in chapters 5-7 use the coordinates u and v for on-panel coordinates. They are based on the 
gauss points on a particular panel. Thus, they don't change with our shape parameters. 


^. = — = 0 (3. 21 ) 

dri dri 


77 , £ - Transformed wing coordinates The square rj, £ (-1.1) planform doesn t 
change with wing geometry changes. 
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_ o 

dri dri 


(3.22) 


- True Planform Coordinates These do change with wing geometry. 


y = Ul + 1) = + D 


(3.23) 


where 6 is the single-wing span, 5 is the single-wing area, and AR is the aspect ratio 
of the single- wing/ 


dy_ _ (1 + n) / AR 
dS ■ 4 V 5 


(3.24) 


dy (l+v) [J_ 
dAR 4 V Ail 


_ o 

dA dA 


(3.26) 


where A is the wing taper ratio. 

x is a function of the local chord length, so and d ^ ] are derived first. 


C r = 




6(1 + A) 1 + A 


(3.27 


dC r _ 1 

dS - 6(1 + A) 


dC r -1 1 / S 

dAR ~ (1 + A) AR V AR 
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^ = 0 (3.30) 

d A 


dCV -2 / 5 

dA “ (1 + A) 2 V Ai? 

(3.31) 

C(y) = C r l-i(l-A)(l + n) 

(3.32) 

where C(y) is the local chord length as a function of the span- wise 

coordinate, y. 

dC(y) dC r C(y ) 

dS dS C r 

(3.33) 

dC{y) dC r C(y) 
dAR ~ dAR C r 

(3.34) 

dC(y) dC r C(y) _ Q 
d\ ~ d\ C r 

(3.35) 

dC(y) dC r C(y) C r * 

d\ dX C r 2 

(3.36) 


Now, back to x: 


x = 



^ C(y) + y tan A 


dx 

dS 




dC{y ) 
dS 


dy . 
— tan A 
dS 


dx 

dAR 



n dew 

4 7 dAi? dAi? 


(3.37) 


(3.38) 


(3.39) 


S.O S*n»titivity Equation* 


27 


(3.40) 


dx _ y 
d\ cos 2 A 


J£ = fl±l. A^M + ^tanA 

d A V 2 


4 y cf A + d\ 


(3.41) 


w 


., Wy - Gauss weights. These don't change with respect to geometry. 


r,(x,y) - Generalized mode shapes. As given in equation 2.5. 


7i(*> y) 


* \ m f y Y 


L £ max / V 2 /max 


(3.421 


where m = 0 to mw and n = 2 to nw + 2. N is looped over first. Thus, 


~n{2,y) = 


X \°{ y Y 


Umax 


3.43) 


72(x,y) 


x \V y A 3 


Umax , 


(3.44 ) 


If m = 0, 


dy _ ny n 1 dy 
dr i V max 


(3.45) 


otherwise ( m / 0 ), 


dy 

dri 


xm-iyn-i 

y m j iti 
x maximal 


dx dy 

my- b nx — 

dr i dr i 


(3.461 


,/j - First Jacobian 


. The Jacobian of the x, y to coordinate transformation is 
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J i 


5 

4 



^(1 + A)J 


Then, 


dJi _ 1 / (1 ~ A) \ _ S (l ~ A ) 

~dS ~ 4 V ^(1 + A)/ 4(1 + A )dS 


recalling that — 0 we get 


(3.4S) 


_ h. (3.49) 

dS 5 


dJi 

d.4J? 


= 0 


(3.50) 


^1=0 (3.51) 

d A 

jji Sr l (3.52) 

dA 2(1 + A) 2 


J 2 - Second Jacobian. The panel based interpolation schemes discussed in chap- 
ters 5-7 use rectangular sub-sections of the square 17 , £ wing defined such that the local 
corrdinates (u. v) = (-1,1) The Jacobian of the rj to u, u coordinate transformation is 


= (3.53) 

S’u, 


where is the panel area measured in the ( 17 , 0 system and 5 u ,u is the panel area 
measured in the (u,u) system which alwavs equals 4. 

The derivatives of J 2 are 
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3-4 Deflection Sensitivities 


These results can be used to calculate a variety of other useful results. For instance, 
equation (2.4) can be differentiated and combined with the equations in section 3.4 to 
produce the sensitivity of the deflection of any point of the wing to changes in the shape 
parameters. Recalling 


np 

W{x.y) = ^7,(x.y)C, 
1=1 


(3.55) 


with 

7i( x -y) = ~ ) m (~ ) n < 3 56 ' 

£max Umax 

where m varies from 0 to mw and n varies from 2 to nw + 2 with n varying more rapidly 
than m. Then 


dW(x,y) 

dr i 


mw nw J f2 


-E 


mc l 

X 

m — 1 

y 

£ max 

_ ^max J 


_ Umax M 


n dx_ 
dr i 


+ 


\dci_ 

X 

m 

y 

n 

1 dr i 

max 


_ y mar _ 




m 


n-1 

nc l 

X 


y 

dy 

ymaz 

_ % max 


_ Umax 

dri 


(3.5 


i ) 
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4.0 Chebyshev Pressure Representation 


The wing pressure field is represented first by a tensor product of Chebyshev polyno- 
mials. They were chosen for their ability to accurately fit a curve with a small number of 
terms, due to their orthogonality properties. The pressure distribution can thus be written 

<3 P 

= <41) 

g=0 p = 0 

Here T p is the Chebyshev polynomial of order p. 

T p (x) = cos(ncos~ l (x)) (4.2) 

Thus, 

T 0 (x) = 1. (4-3) 

Ti(x) = x. (4.4) 

T-i(x) = 2x 2 — 1, (4.5) 

and in general 

T n+1 =2xT n (x)-T n _ 1 (x). (4.6) 

A quick comparison of Eqns. 4.1 and 2.9 shows that the generic 3 J is 

^■(17.0 =r,(i7)T f (0 (4.7) 

4-1 Chebyshev Results 

The aeroelastic code does an excellent job of calculating the converged wing loading 
and deflections for a particular flight condition. The sensitivity derivative calculations do 
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an excellent job of predicting shape sensitivities. Slight difficulty has been encountered in 
the exact calculation of derivatives with very small values. It was modeled structural!} 
using a box beam model detailed in Ref. 10. The material properties were E n = E 21 = 

6.S9 x 10 10 , G 1 2 = 2.65 x 10 l °, and u l2 = 0.3. 

The behavior of the chebyshev aerodynamic coeffients, a, are of interest because 
several perturbed runs of the program are needed to produce finite difference derivati\es 
for later use by the sensitivity program. As Figures 8-9 show, the overall trend of these 
aerodynamic coefficients is smooth, but with significant sharp “wiggle" in these curves. 
This wiggle makes accurate finite differencing problematical on the affected a s. These 
curves are the worst cases; many of the other curves are much smoother. 

To overcome this problem, a higher-order finite difference scheme using a large step 

size was used. This scheme 

/ '( x ) « —[/{x - 2Ax) - Sf(x - Ax) + 0 f(x) + 8/(x + Ax) - f(x + 2Ax)] (4.8) 

12 h 

has an accuracy of 0(Ax 4 ). 

Note, these figures are for single elements of the sixty element {a} vector. These 
graphs can be used to estimate the best step size for a single element, but not for the 
entire vector. In an attempt to discover an optimal step size to maximize the accuracy of 
the vector as a whole, the L-2 norm of the difference between two finite differences was 
minimized. In this case, for each of a wide variety of step sizes, a forward and central first 
derivative finite difference was calculated with the central difference used as a reference. 
The logarithm of the average of the square of the difference of these two values was plotted. 
The step size that produced the minimum difference was chosen as our optimal step size 
for use in eqn. 4.7. Note, this is not in any way a percentage error, but instead the average 
absolute squared error. Figure 10 shows the variation of this pseudo-error with step size 
in a and aspect ratio. Similar curves were generated for the other three n s. 
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Figure 8. Pressure Expansion Coefficient a(l,l) vs. Sweep 
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Figure 9. Pressure Expansion Coefficient a(3,4) vs. Sweep 
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Figure 10. Finite Difference Pseudo-Errors 
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The individual deflection coefficients, C. show some wiggle, similar to the aerodynamic 
terms, but the effect is much smaller. Since these derivatives are only used to confirm the 
sensitivity results, a step size study was less important.. Nevertheless, one was performed 
to insure that these derivatives would be as accurate as possible. 

Derivatives of the stiffness matrix are also calculated by finite difference for use in the 
sensitivitv code. A step size study identical to the aerodynamic code study was undertaken. 
While the term “L-2 Norm” is not strictly accurate for a matrix operation, the process 
was the same. The difference between each term of a central and a forward differenced 
derivative was squared and summed. The behavior of this pseudo-error measurement is 
also plotted in Figure 10. The optimal step size was used by a central difference in the 
actual sensitivity calculation. 

The variation of the trim angle of attack with respect to the wing area is shown 
in Figure 11. The solid line shows the converged results from the iterative aerodynamic 
and structures combination. The various dashed lines show the variation predicted by the 
sensitivity derivatives at the different base configurations. The prediction goes through 
the converged value at the base geometry and is linear with a slope equal to the sensitivity 
derivative. The desired result is for this line to be tangent to the converged data curve. 

Similarly, the sensitivity of the trim angle of attack to changes in the wing aspect ratio 
is shown in Figure 12. The solid line shows the converged iterative results and the dashed 
lines show the predicted variation by having a slope equal to the calculated sensitivity 
derivative. 

Figures 13 and 14 show the converged and predicted values for the angle of attack 
variation with respect to taper ratio and sweep. It is obvious from Fig. 7 that the obtained 
value of the sensitivity of the angle of attack with respect to the taper ratio is not very 
accurate corresponding to the taper ratio values of 0.5 and 0.8. However, note that the 
value of the converged angle of attack is almost insensitive to the variation in taper ratio 
at those values of taper ratio. The inaccuracy in the present results can be attributed to 
the numerical problems associated with determining derivatives that are almost zero. 

Also of interest are deflection sensitivities. Figures 15-18 show the predicted and 
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Figure 11. Trim Angle of Attack vs. Wing Area, Chebyshev Pressure Representation 
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Figure 12. Trim Angle of Attack vs. Wing Aspect Ratio, Chebyshev Pressure Representation 
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Figure 13. Trim Angle of Attack vs. Wing Sweep, Chebyshev Pressure Representation 
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Figure 14. Trim Angle of Attack vs. Wing Taper Ratio, Chebyshev Pressure Representation 
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Figure 15. Deflection Coefficient Cl vs. Wing Area, Chebyshev Pressure Representation 
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Figure 16. Deflection Coefficient Cl vs. Wing Aspect Ratio, Chebyshev Pressure Representation 
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Figure 17. Deflection Coefficient Cl vs. Wing Sweep, Chebyshev Pressure Representation 
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Figure 18. Deflection Coefficient Cl vs. Wing Taper Ratio, Chebyshev Pressure Representation 
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actual variation of a few of the deflection expansion coefficients. Figure 15 is the worst 
of all of the plots in this report. Note the relative scaling of the axes in this plot. The 
sensitivity is very near zero. Perhaps more important is the sensitivity of the deflections 
themselves. Figures 19-22 show the sensitivity of the leading edge tip deflection (measured 
in meters). The sensitivity of the coefficient C\ with respect to wing area is not very 
accurate (see Fig. 15) as C\ is almost insensitive to the variation in the wing area. The 
inaccuracies obtained in calculating this sensitivity do not however affect the sensitivity 
of the tip deflection with respect to wing area. This can be observed in Fig. 17. From 
Figs. 17-20 it is clear that the present formulation yields very accurate shape sensitivities 
for the aeroelastic tip deflections and can be used in optimization studies. The coefficient 
errors largely disappear in the integrated quantities. 


Table 1: 

Comparison of Finite Difference and Analyt: 

ic Logarithmic Derivatives 


Term 

Finite Difference 

Analytic 

% Error 

da 

dS 

da 

dAR 

- 1.188194 

-1.129775 

4.917 

-0.366320 

-0.463287 

29.906 

da 
d \ 

-1.85613xl0 -2 

-2.40882xl0 -2 

29.777 

da 

d\ 

0.098611 

0.126455 

28.236 

dc\ 

dS 

-5.83136xl0" 4 

-9.77324xl0 -3 

40.333 

dc i 

dAR 

0.122537 

0.123591 

0.8581 

dci 

dX 

0.332564 

0.329906 

-3.957 

dc\ 

d\ 

0.122932 

0.112007 

-8. 887 

dTipDtf 

dS 

0.984423 

0.978191 

0.2812 

dTipDt f 
dAR 

1.888806 

1.916328 

1.457 

dT ipDe f 
d\ 

0.273598 

0.265756 

2.866 

dTipDef 

0.262560 

0.258709 

1.467 


a a 

These errors are largely numerical in origin. Variables with very small logarithmic 
derivatives will be difficult to differentiate numerically regardless of the scheme used 21 . 
Table 1 shows a variety of variables, their logarithmic derivatives, and the error in calcu- 
lating them. As can be seen from the table, some of the intermediate results with very low 
logarithmic derivatives show as much as a 40% error. However, sensitivities of the final 
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Figure 19. Tip Deflection vs. Wing Area, Chebyshev Pressure Representation 
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Figure 21. Tip Deflection vs. Wing Sweep, Chebyshev Pressure Representation 
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Figure 22. Tip Deflection vs. Wing Taper Ratio, Chebyshev Pressure Representation 
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integrated tip deflection results show a maximum of a 3% error. 

The results show good accuracy for integrated quantities such as tip displacements but 
less accuracy for individual displacement coefficients or trim angle of attack. In general, 
the accuracy decreases noticably when the size of the derivative decreases. 

The global sensitivities of the aerodynamic coeffients were also studied. Equation 3.4 
was used to calculate these sensitivities. As with the deflection coefficient sensitivities, 
it was difficult to accurately model the sensitivities of individual coefficients. A wide 
variety of coefficients were examined for each of the four shape parameters. Figures 23-30 
are representative of the overall results. Figures 23-24 show the mediocre modeling of 
da/dS. The similarly bad modeling of da/dAR is shown in figures 25-26. Despite more 
eratic behavior in the variation of the a’s with respect to sweep and taper, the modelings of 
da/dA and da/dX are much better. These are shown in figures 27-2S and 29-30 respectively. 
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Figure 23. a(l) vs. Wing Area, Chebyshev Pressure Representat.on 
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Figure 25. a(2) vs. Aspect Ratio, Chebyshev Pressure Representation 
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Figure 29. a(15) vs. Taper Ratio, Chebyshev Pressure Representation 
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5.0 Constant Local Pressure Representation 


The globed Chebyshev pressure representation will not adequately capture small, lo- 
cal, details in the pressure field unless a very large number of terms are used. Thus, 
an edternative scheme was employed. The pressure is represented by using interpolation 
polynomials that have local support over the wing panel. It is expected that a larger 
aerodynamic coefficient vector, {a}, will be necessary for similarly accurate results. This 
relationship is examined. 

Any number of polynomials are available for this purpose. To begin with a scheme 
that represented panels each having constant pressure was used. Given a square panel 
defined as (u,v) = (-1, 1), and the values at the four corners of the panel as 6 0 o- &oi- &io- 
b n . Then the interpolated value p at any point (u. v) on that panel is given by: 

p(u.t') = ^[l 1] 

This can be rearranged by multiplication to become: 


&oo 

&io 


froi 

b n 


(o.i : 


p(u,v) = - [1 1 1 1] 


b oo 

boi 

&io 

l 6n / 


> 


( 5 . 2 ) 


In other words, the pressure on a panel is represented by the average of the pressures 
at the four corners of that panel. The symbol R,(u,v ) is used in the following chapters 
to represent the interpolation polynomials used by the bilinear and biquadratic interpo- 
lation. For completeness sake, the constant pressure panel representation s interpolation 
polynomials are 
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R t {u.v) = i 


(5.3) 


where i varies from 1 to 4. 

This form is suitable for the integration of the forces. If we call the right hand side 
our aerodynamic coefficients (actually true pressures), {a}, we get 

p{u,v) = {i?} r {a} (5.4) 

This pressure representation can be used in eqns. 2.10 and 2 . 12 . These equations 
allow us to produce the [.4], aerodynamic kernal matrix, for the constant pressure panel 
case. A typical term in [Ac], our constant panel aerodynamic matrix, is: 


i A c}= f f ^( u ' l ’)7.(2:-y)l>/i('7,Om2(u.i')|du<ir 
». J — i J 


(5.5) 


panels 


where «/ 2 ( u , t; ) is the Jacobian of the transformation from the (/?,£) to (u, v) coordinate sys- 
tem, and the subscript C indicates the constant pressure representation of the aerodynamic 
pressure over the wing panels. 

[Ac] has the following properties: 

1 ) Each row corresponds to a particular Q t (corresponding to 7 ^). 

2 ) Each column corresponds to the contribution of a particular point pressure. 

3) Each wing panel will add to 4 columns in [.4], For instance, panel will add to 
each element u z” in the column corresponding to the upper left node: 



71(^0 ki(^Oll^2(u.u)l du dv 


(5.6) 


The constant-panel version of {£} is very similar. Each panel contributes to the four 
elements of the vector corresponding to its four corner nodes. 
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(5.7) 


L c = Y, f f R{u ' v) du dv 

panels 1 ^ 

A Gauss integration scheme was used to perform the double integrals in the [A c ] and 
{Lc} equations. 

The first step was to verify [Ac] and {Lc}- Recalling that Q, = [A c ]{«} and 
Lift = {Lc) T { a ) gives us an easy way to test the terms. Comparing Q, and the lift 
to results earlier obtained from using global Chebyshev polynomials eventually gave very 
good results. 

An evenly distributed rectangular grid was first used. A computer graphics visualiza- 
tion of the error field confirmed that significant error was concentrated in the region very 

near the leading edge. 

A revised gridding scheme was devised. This had two refinements. The inability of 
FAST to produce pressure data at the wing edges was initially a problem. Lnder the 
original local bilinear scheme, panels were assembled away from the edge, resulting in a 
100% error band around the boundary of the wing. The current scheme puts the panels 
nearer to the wing edges and also specifies a zero pressure along the edges, allowing panels 
to be placed along the wing edges. 

Second, a variable spacing technique was used to put improved resolution where it 
was needed. This used a modified half-Gauss-Lobatto grid. The phase and frequency of 
the cosine term in the Gauss-Lobatto equation was changed to give high resolution at the 
wing leading edge and tip, and low resolution along the trailing edge and wing root. 


, 7T Z — 1 

£(*') = -2cw(- — ) - 1 


1 5.8) 


n (j) = -2cos(- J ~^— A - i ~ 1 


).0l 
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This scheme captures the leading edge peak of the pressure field much better than 
the earlier scheme employing a uniformly distributed grid and also better than the global 
chebyshev scheme employed in Chapter 4. 

Then, the derivatives of [A] and {£}, taking into account all the zero coordinate 
derivatives, are: 


dAc 

dri 


£ 5 > ' Z WyJ 2 
kgx kgy 


( 


d^f » _ d 3 1 

dF l J ' +1 ’Tr 7 


(5.10) 


and 


dL 

dr i 


kgx kgy 


(5.11) 


where kgx and kgy are the number of gauss points in the x and y directions. (Typically, 
kgx -- kgy = 2). Note that dJo/dri = 0. 

Tables 4 - 13 (in Chapter S) show the results from the constant pressure panel versions 
of the analysis. Comparisons with the other approaches show very close results. 
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6.0 Bilinear Local Pressure Representation 


Next, a bilinear scheme was used to represent the wing pressure. Given the same 
square panel defined as {u, v) = (-1,1), and the values at the four corners of the panel as 
b 00 , &oi, 610, &11. Then the interpolated value p at any point (u,u) is given by: 

p(u,v) = i[l - u 1 + u] 

This can be rearranged by multiplication to become: 


boo 

^10 


&01 

bn 


1;: 


(6.1: 


p{u. v) = - [(1 - u)(l — v) (1 — u)(l+u) (l + u)(l — v) (1 + u)(l + i>)] 



( 6 . 2 ) 


iL 


If we call the left hand vector R. our interpolation vector, and the right hand side is 
our aerodynamic coefficients (actually true pressures), {a}, we get 

p(u, v) = {i?} 7 '{a} ( 6 . 3 ) 

where a typical term in [As z,]. our bilinear aerodynamic matrix, is: 

[a B l] = 51 / / R ( u ^ v hi( x 'y)\ J ^ r i'^\ J ^ u ' v ^ dudv ( 6 - 4 ) 

panels ^ ^ 

where / 2 (u, v) is the Jacobian of the transformation from the (r?,£) to (u,u) coordinate 
system, and the subscript BL indicates the bilinear representation of the aerodynamic 
pressure over the wing panels. 
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[.4s/,] has the same properties discussed for [-4c]- For instance, panel "j will add to 
each element "i" in the column corresponding to the upper left node: 

/:/: i(l - u)(l - v) 7 ,{r).0 \Ji{r),0\\Hu>v)\ du dv ( 6 - 5 ) 

The bilinear version of {L} is very similar. Each panel contributes to the four elements 
of the vector corresponding to its four corner nodes. 

Lql = E fj\ R(u,v) |Ji(j 7 ,OII^ 2 (u,u)| du dv (6.6) 

panels ^ 

A Gauss integration scheme was used to perform the double integrals in the Abl and 
Lbl equations. 

The next thing to test was the derivatives of [As/,] and {£}. This was accomplished 
by comparing the analytic derivatives to finite difference derivatives calculated by recalling 
the subroutines with slight changes to the shape parameters. 

The new bilinear-interpolation versions of the terms {d{a} / da} , [d{a }/d{C}}. and 
{d{a}/dr ( } remained. Each was calculated by appropriate finite differences. Then, the 
sensitivity code had to be modified to allow easy switching between the two (and other 
future) schemes. 

R(u,v) - Bilinear interpolator. R has 4 different values, one at a time, depending 
on which term of [A] is under consideration. 

R l = 0.25(1 - u)(l - u) = 0.25(1 - u - v + uv) (6.7) 

R 2 = 0.25(1 + u)( 1 - v) = 0.25(1 + u-v-uv) (68) 

= 0.25(1 - u )( 1 + v) = 0.25(1 - u + b- uu) (6.9) 
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R 4 = 0.25(1 +u)(l + v) 


( 6 . 10 ) 


Since the derivatives du/dr ( and dv/dr( are zero, the derivatives of R, is also zero. 


6.1 Bilinear Pressure Interpolation 


The motivating factor for attempted local or panel based schemes was to try to 
improve the accuracy of the local derivatives of the aerodynamic loads. A major factor in 
the inaccuracy of the Chebyshev approach’s aerodynamic derivatives were the "wiggles" 
in the data being differentiated. 

A study was undertaken to determine if the new aerodynamic coefficients, actually 
point pressures, were any smoother in their behaviors. Each of the four shape parameters 
was varied and pressures at three points on the wing were stored. These points were the 
nodes nearest to the wing leading-edge- root, the wing center, and the wing trailing-edge- 
tip. Nearly ail of the data was if not linear, at least very smoothly varying. Figures 31 - 34 
show the four worst cases. Only the trailing-edge-root pressure vs. Sweep plot indicates 
a potential problem. The other sweep plots were perfectly smooth. The pressures shown 
are measured in Pascals. 

6.2 Bilinear Results 

The variation of the trim angle of attack with respect to the wing area is shown 
in Figure 35. The solid line shows the converged results from the iterative aerodynamic 
and structures combination. The various dashed lines show the variation predicted by the 
sensitivity derivatives at the different base configurations. The prediction goes through 
the converged value at the base geometry and is linear with a slope equal to the sensitivity 
derivative. The desired result is for this line to be tangent to the converged data curve. 
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Figure 33. Trailing Edge-Tip Pressure vs. Sweep 
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Figure 35. Trim Angle of Attack vs. Wing Area, Bilinear Pressure Representation 
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Similarly, the sensitivity of the trim angle of attack to changes in the wing aspect ratio 
is shown in Figure 36. The solid line shows the converged iterative results and the dashed 
lines show the predicted variation by having a slope equal to the calculated sensitivity 
derivative. 

Figures 37 and 38 show the converged and predicted values for the angle of attack 
variation with respect to taper ratio and sweep. It is obvious from Fig. 38 that the 
obtained value of the sensitivity of the angle of attack with respect, to the taper ratio is 
not as accurate as one would like. However, note that the value of the converged angle of 
attack is almost insensitive to the variation in taper ratio at those values of taper ratio. The 
inaccuracy in the present results can be attributed to the numerical problems associated 
with determining derivatives that axe almost zero. 

Also of interest are deflection sensitivities. Figures 39-42 show the sensitivity of 
the leading edge tip deflection. From these figures it is clear that the present formulation 
yields very accurate shape sensitivities for the aeroelastic tip deflections and can be used in 
optimization studies. The coefficient errors largely disappear in the integrated quantities. 

These few errors are largely numerical in origin. Variables with very small logarith- 
mic derivatives will be difficult to differentiate numerically regardless of the scheme used 
(Haftka 23 ). The final comparisions between the current versions of both schemes is also of 
interest. The following chart. Table 2, tompares the various results. The finite difference 
derivatives are all forward derivatives where the step size was decreased until the derivative 
converged: 
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Figure 36. Trim Angle of Attack vs. Wing Aspect Ratio. Bilinear Pressure Representation 
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Table 2: Comparison of Bilinear and Chebyshev approach Derivatives 


W 


Term 

Bilinear 

Analytic 

Chebyshev 

Analytic 

Bilinear 
Fin. Dif. 

Chebyshev 
Fin. Dif. 

trim 

6.2686 

6.0987 



TD m 

0.4086 

0.3879 



TDiv/a 

0.6179 

0.5915 



da/dS 

-0.710 

-0.689 

-0.6S6 

-0.725 

dTD/dS 

0.040 

0.038 

0.040 

0.038 

da/dAR 

-0.766 

-0.775 

-0.579 

-0.596 

dTD/dAR 

0.209 

0.274 

0.209 

0.195 

da j dA 

-0.056 

-0.051 

-0.032 

-0.040 

dTD/d A 

-7.2xl0 -3 

-6.7xl0" 3 

-7.2xl0 -3 

-6.SxlO -3 

da/d A 

-0.359 

-0.294 

-0.188 

-0.226 

dTD/d X 

0.213 

0.266 

0.227 

0.212 


* Tip Deflection 
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Figure 38. Trim Angle of Attack vs. Wing Taper Ratio, Bilinear Pressure Representation 
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Figure 40. Tip Deflection vs. Wing Aspect Ratio, Bilinear Pressure Representation 
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Figure 41 . Tip Deflection vs. Wing Sweep, Bilinear Pressure Representation 


6.0 Bilintar Lock] Pressure R«pr**« ntlt * on 



Tip Deflection 



Wing Tip Deflection vs. Taper Ratio 
Bilinear Pressure Representation 


Figure 42. Tip Deflection vs. Wing Taper Ratio, Bilinear Pressure Representation 
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7.0 Biquadratic Local Pressure Representation 


A biquadratic pressure panel representation was developed next. Again, a panel was 
designated as having (tt.r) = (-1. 1). For this scheme, the panel has nine input pressure 
nodes. 

The interpolation polynomials R , were derived from Langrangian polynomials. The 
two dimensional R's are products of the one dimensional quadratic Lagrangian polynomi- 
als. These Lagrangian polynomials with an arbitrarily positioned middle node are: 


Li(u) 


1 (u — t/2)(n — 1) 

2 (1 + U 2 ) 


(7.1) 


_ (u + 1 Xu - 1 ) 

(ti2 + l)( u 2 — 1) 


(7.2) 


__ 1 ( U - U 2 )(» + 1 ) 

= 2 (l-u 2 ) 


(7.3) 


where 

v is the location to be interpolated to and 
u 2 is the coordinate of the middle node. 

Next, we define our two dimensional interpolation polynomials as follows: 


Ri(u,v) = Li(u)Li(v) 

(7.4) 

R 2 {u,v) = I 2 (w)Ii(t>) 

(7.5) 

H 3 (u,i>) = L 3 (u)Li(it) 

(7.6) 
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.R 4 (u, v) = Li(u)L 2 (v ) 


(7.7) 


Rs(u,v) = L 2 (u)L 2 {v) 

(7.8) 

Rg(u.v) = L 3 (u)L 2 (v) 

(7.9) 

Ri(u,v) = Li(u)Li(v) 

(7.10) 

Rs(u,v) = L 2 (u)L 3 (v) 

(7.11) 

Rg(u,v) = L 3 (u)L 3 (v) 

(7.12) 


This development allows us to use the same half-Gauss-Lobatto grid scheme developed 
for the constant and bilinear panels. Initially, this didn’t work properly. The biquadratic 
interpolation calculated a significantly higher lift for a given deflection, leading to a sig- 
nificantly lower required trim angle of attack, roughly half of what was predicted by the 
other methods. 

Some graphical visualization revealed the problem. The pressure field calculated by 
FAST is shown in figure 33. The field is heavily dominated by a huge spike along the 
leading edge. The constant and bilinear fits were able to deal with what is effectively an 
impulse in the field. The biquadratic representation could not. As it tried to fit a curve 
through the zero pressure on the leading edge, through the huge value at the next node, 
and back to a fraction of that value on the third node, it encountered numerical difficulties. 

The solution to the problem was to reduce the gradient of the pressure field at the 
leading edge. For this purpose, the edge values were set equal to the values at the near- 
est internal node. Thus, the edge pressure gradient was set to zero. The biquadratic 
interpolation then worked correctly. 

To understand the magnitude of the net change in the pressure field, the bilinear and 
constant convergence tests were rerun. The converged value of trim angle of attack 
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Figure 43. Wing Pressure Field 
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I 

m 

was reduced by about 2%. The change is small because the area of the panels effected is 
extremely tiny. 

7.1 Biquadratic Results 

Figures 44-47 show the biquadratic pressure representation results for trim angle of 
attack and the sensitivity derivatives. Figures 48 - 51 show the tip deflection values and 
their sensitivity derivatives. These plots are nearly identical to those produced by the 
other pressure schemes. 

Chapter 8 discusses the relative performance of the schemes in detail. 

Equation 3.4 was also used to try to represent the global sensitivities of the point 
pressures to the shape parameters. It was anticipated that as with deflection coefficients 
it would be difficult to model the individual sensitivities of the 1600 coefficients. As with 
the Chebyshev version of da/dr/, the results were mixed. Figures 52-59 are representative 
of the large number of point pressures examined. 

Figures 52-53 show da/dS at two points along the leading edge. Figures 54-55 show 
da/dAR's mediocre performance. As with the chebyshev approach, da/dA and da/dX per- 
form much better. Representative plots are shown in figures 56-57 and 58-59 respectively. 


• ft 


T.O Biquadratic Local Prua«r« Representation 



i 



Figure 44. Trim Angle of Attack vs. Wing Area, Biquadratic Pressure Representation 
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Figure 45. Trim Angle of Attack vs. Wing Aspect Ratio, Biquadratic Pressure Representation 
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Figure 46. Trim Angle of Attack vs. Wing Sweep, Biquadratic Pressure Representation 
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Figure 47. Trim Angle of Attack vs. Wing Taper Ratio, Biquadratic Pressure Representation 
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Figure 48. Tip Deflection vs. Wing Area, Biquadratic Pressure Representation 
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Figure 49. Tip Deflection vs. Wing Aspect Ratio, Biquadratic Pressure Representation 
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Figure 50. Tip Deflection vs. Wing Sweep, Biquadratic Pressure Representation 
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Figure 51. Tip Deflection vs. Wing Taper Ratio, Biquadratic Pressure Representation 
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Figure 53. L.E. Tip Pressure vs. Area, Biquadratic Pressure Representation 


^ ® Biquadratic Local Pressure Representation 


03 



Leading Edge Root Pressure 



Figure 54. L.E. Root Pressure vs. Aspect Ratio, Biquadratic Pressure Representation 
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Figure 58. L.E. Root Pressure vs Taper Ratio, Biquadratic Pressure Representation 


7.0 Biquadratic Local Pressure Representation 





8.0 Method Comparison 


In this research, a variation of Sobieski s Global Sensitivity Equations is implemented 
to obtain the global sensitivity of the static aeroelastic responses. The scheme is indepen- 
dent of the analysis code used to obtain aerodynamic data. 


8.1 CPU Time Comparison 


Each method produced very similar final results. The remaining question is how well 
they performed costwise. The following chart shows a comparison of CPU times on the 
VAX mainframe for a complete analysis run. Here, Nx is the number of pressure nodes in 
the x, or chord-wise, direction. Ny is the number of pressure nodes in the y, or span- wise, 
direction. 


Table 3: CPU Times 


Case 

HR 

CPUTime 
: MIN : SEC 

Chebyshev 

1 : 31 : 55.31 

Constant Nx=Ny=15 

1:11: 42.02 

Constant Nx=Ny=39 

1 : 24 : 41.62 

Bilinear Nx=Ny=15 

1 : 20 : 52.99 

Bilinear Nx=Ny=39 

1 : 23 : 00.35 

Biquadratic Nx=Ny=15 

1 : 14 : 17.58 

Biquadratic Nx=Ny=39 

1 : 24 : 41.40 
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Thus, the differences in CPU run times for the various schemes and resolutions vary 
by less than 20%. 


8.2 Performance Comparison 


The next issue for the local schemes is just how much resolution is necessary. This is 
not a simple issue. 

Tables 4-33 show' the various converged results for all three local pressure schemes. 
A wide range of resolutions is investigated. The tables also list each value’s percentage 
difference from the highest resolution grid’s value. These percentage differences give an 
indication of how well a particular grid has done at producing a converged result. The 
sensti\it\ derivative charts also give the logarithmic derivative of the highest resolution 
grid's value. This value gives some indication of how accurately it is possible to model the 
derivative. 

Some grid conclusions are valid across all of the cases. It is obvious that _V r has a 
much more significant effect on the solution than A y . This is because the main variation in 
the pressure field is in the x. or chordwise, direction. A certain minimum resolution in the 
y. or span wise, direction seems necessary, perhaps N y = 15. After that, increasing spanwise 
resolution accomplishes little compared to a similar increase in chordwise resolution. 

Since the difference in CPL times between high and low resolution cases is not very 
large, in most cases running the program at the highest resolution available makes sense. 
If how’ever, one wishes to minimize the CPU time used, then some choices must be made. 
One has to decide how r close to a converged value one is willing to be. The convergence 
tables give some qualitative feel for what grid spacings are necessary to be as close as you 
desire to convergence for a particular technique. 

For instance, assume a 10 % difference from the highest resolution value is permissable 
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Li 


for a certain application. A look through the Biquadratic results shows that a 15x15 grid 
will provide this accuracy. The bilinear approach requires a 15x5 grid to get this result. 
The constant pressure panel approach requires a 15x15 grid. 

The choice of which scheme to use remains. All four schemes work well. Of the three 
local schemes, the higher order schemes converge faster than the lower schemes, but the 

difference isn't overwhelming. Since the CPU times for the three schemes are comparable, 
the biquadratic scheme produces the best results for similar computational effort. However, 
the polynomial overshoot problems encountered in this approach could be problematical 
for other cases. The bilinear approach is more robust and only slightly slower to converge. 


Trim Angle of Attack, Degrees 

Local Constant Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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Table 4. Trim Angle of Attack, Local Constant Pressure Representation 
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Leading Edge Tip Deflection 

Local Constant Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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da/dS 

Local Constant Pressure Representation with Zero Gradient Edge: 
(percentage difference from 39x39 value) 
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Table 6. da/dS , Local Constant Pressure Representation 
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Logarithmic Derivative: 0.87825 


da/dAR 

Local Constant Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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Table 8. da/dAR , Local Constant Pressure Representation 
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(3.90%) (2.94%) (9.00%) (6.93%) (0.0%) 

Logarithmic Derivative: -0.37747 
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Local Constant Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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Table 10. da/dA, Local Constant Pressure Representation 
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Table 11. dTip Deflection, 'd A, Local Constant Pressure Representation 


$.0 Method Compiriion 


Logarithmic Derivative: 0.20771 


da/cR 

Local Constant Pressure Representation with Zero Gradient Edg< 
(percentage difference from 39x39 value) 
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Table 12. da/dX , Local Constant Pressure Representation 


B.O Method Companion 


Logarithmic Derivative: -2.49x10 
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Table 13. dTip Deflection/d A, Local Constant Pressure Representation 


8.0 Method Comparison 


Trim Angle of Attack, Degrees 

Local Bilinear Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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Table 14. Trim Angle of Attack, Local Bilinear Pressure Representation 


ft.O Method Comparison 
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Table 15. Leading Edge Tip Deflection, Local Bilinear Pressure Representation 


• .0 Method Comparison 
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(1.63%) (0.202%) (0.052%) (0.008%) (0.0%) 
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Local Bilinear Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 


ON 

CO 

II 


2 


NO 

, s 

00 ^ 


^ P' 

CO 


O >© 

CO 

^ o 5 

ON 

<N 


CO g' 
LO 

ON rH 

LO CM 

CM 

NO LO 

NO 

ON 

NO P 

NO LO 

NOo 

d 

i 


d 

* • 
90 

9 S 


o 

00 

CM 

NO 

NO 

o 


V© 

o 


LT) 

co 

II 

Sn 


2 


r — 1 


NO 

x v 

in 

o 

in 

. o 

NO 

•^p 

00 

v o 

04 

S.O 

in 

0 s 

04 


04 


ON 


CO 


ON 


04 

o 

CO 

ON 

04 

00 

LO 


LO 

LO 

NO 

00 

NO 

NO 

NO 

ON 

NO 


NO 


NO 

o 

NO 

O 

d 

i 


d 

rH 

d 

i 

O 

d 

i 

d^ 

d 

i 




CM ^ 

H vO 

o g; 
co £? 
nD 

d p 


CM 

CO 

LO 


vO 

o'" 


d 


CM 3 * 

ON o 

LO rH 

NO LO 
• • 


^ P" 

LO 

CM 00 

NO co 

NO o 

• « 

90 


PI cr 

S ©^ 
co co 

NO CO 

NO o 

• • 

P o 


LO 


II 


;n 

2 


NO 



P 1 P" 

co P" 

co P 

co 


LO -V© 

lO o'' 

On 

LO o'' 

o 


CO ©i 

On i 

CM O 

CO rH 

CO 

ON 

LO 2 

LO NO 

NO T-l 

NO rH 

NO 

00 

NO ^ 

Np 

NO o 

NO T - • 

d 

i 


d d. 

d 6 


?o 


LO 

II 


2 


ON 


CM 


00 

^ D 

LO 

. o 

IN 

N^O 

NO 


ON 


CM 



cO 

00 

O'' 

(N 


LO 


CM 

00 

m 

T — 1 

LO 


co 


LO 


NO 

IN. 

NO 

o 

NO 

NO 

NO 

uh> 

NO 

O 

NO 

O 

NO 


NO 


d 

i 


d 

i 

r-H 

d 

d^ 

o 

i 

d^ 

d 

i 

d 


LO 

II 

X 


2 



Table 16 . da/dS, Local Bilinear Pressure Representation 


ft.O Method Comptriion 


115 


Logarithmic Derivative: -1.12328 
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Table 17. dTipDeflection / dS, Local Bilinear Pressure Representation 


1.0 M • t h od Comparison 
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ogarithmic Derivative: 0.87827 


lx ^ 

LO \C 
Tb £ 
LO P 

o ^ 

T— ( xO 

vO 

00 p 

CO sC 
IX <£ 

00 P 

o ^ 

LO s° 
Ov 

o ^ 

fcc? 

LO Z* 

LO o 

LO VO 

LO VO 

VO rH 

d rL 

i 

o rL 
1 — 

o rL 

i — 

d CS 

dB 


rtf 

<N £ 
VO ON 
LO rH 


O ^ 

r- 1 vP 

VO o' 

<N 

«Op 

o °o 


vO 

VO ^ 
vo I s 
P csi 

o rL 


00 ^ 
00 vO 

o £ 

ID £> 
LO P 
dpv 


X ^ 
CO -c 
LO c 
X ^ 

P 

o rL 


> r- 

5 o 2 

^ XJ '— 

r- » -V fC QJ 

T3 S 
^ c c 
\ <V <D 

X! o£ 


00 ^ 

r— ( vO 

O ^ 
CO »g 
LO ^ 

d lo 


^ ^~v 

\C 

LO ^ 

co s 

LO P 

d vo 


s cr 

ro ^ 

o ^ 

VO .—V 

vO 0 s 
VO o 

rH vO 

vo $: 

vo vp 
o f£ 

LO -p 
rf - 

P IX 

x p 

o P 

^ n 

LO x 

LO 

p ov 

top 

9 a 

o rL 

d rL 

d P 


ix ^ 

\0 
^ cO 

IX p 

P T* 
o »-L 


^ sO 

oo 

co ov 

LO O 

d^ 


C0 x-v 
CO >vO 
Csl O^ 

s® 

LO r— I 


00 ^ 
Tf sO 

00 X 

o LO 

LO^ 

o lo 


(N ^ 

r-H \0 

CM 

LO <N 
LO P 

d ^ 


<N ^ 

ro vp 

vo o^ 

00p 

LO VO 

o rL 

I 


OV 

00 

rH O- 

On P 
LO x 

o rL 


Cs | ^ 
CO sO 
Csl 

OOP 

LO LO 


VO ^ 

rH O- 

t> P 

LO CO 

o rL 


Table 18. da/dAR, Local Bilinear Pressure Representation 


8.0 Method Comparison 


Logarithmic Derivative: -0.56865 
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Table 19. dTipDeflection / dAR, Local Bilinear Pressure Representation 


1.0 Method Comparison 
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da/dA 

Local Bilinear Pressure Representation with Zero Gradient Edges 
(percentage difference from 39x39 value) 
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Table 20. da/dA , Local Bilinear Pressure Representation 


t.O Method Compiriio# 
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.ogarithmic Uenvative: 
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Table 21. dTip Deflection/d A, Local Bilinear Pressure Representation 


• -0 Method Comparison 


Logarithmic Derivative: 0.20764 


da/dX 

Local Bilinear Pressure Representation with Zero Gradient Edgi 
(percentage difference from 39x39 value) 
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Table 22. da/dX , Local Bilinear Pressure Representation 


1.0 Mctkod Compkrii 
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Logarithmic Derivative: -2.49x10" 
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Table 23. dTip Deflection/d A, Local Bilinear Pressure Representation 


1.0 Vfatkod Comparison 
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(1.95%) (0.030%) (0.180%) (0.200%) (0.0%) 

Logarithmic Derivative: 0.27136 


Trim Angle of Attack, Degrees 

Local Biquadratic Pressure Representation 
(percentage difference from 39x39 value) 
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Table 24. Trim Angle of Attack, Local Biquadratic Pressure Representation 


ft.O Method Comptriaot 
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Table 25. Leading Edge Tip Deflection, Local Biquadratic Pressure Representation 


t.O Method Comp* 
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da/dS 

Local Biquadratic Pressure Representation 
(percentage difference from 39x39 value) 
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Table 26 . da/dS , Local Biquadratic Pressure Representation 


• .0 Method Comp&ritoa 
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Logarithmic Derivative: -1.12396 


dTipDeflection/dS 

^ocal Biquadratic Pressure Representatiot 
(percentage difference from 39x39 value) 
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Table 27. dTipDeflection / dS, Local Biquadratic Pressure Representation 


1.0 Method Comptrifoi 


i 


139 


ic Derivative: 0.80755 


da/dAR 

Local Biquadratic Pressure Representation 
(percentage difference from 39x39 value) 
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Table 28 . da/dAR , Local Biquadratic Pressure Representation 


• .0 Matkod Conptruoi 
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dTipDeflection/dAR 

^ocal Biquadratic Pressure Representatior 
(percentage difference from 39x39 value) 
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Table 29, dTipDeflection / dAR, Local Biquadratic Pressure Representation 


§.0 Method Comperiioa 


131 


da/dA 

Local Biquadratic Pressure Representation 
(percentage difference from 39x39 value) 
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Table 30. da/d.\, Local Biquadratic Pressure Representation 


9.0 Method Comparison 
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Logarithmic Derivative: -5.74x10" 


dTipDeflection/dA 

-ocal Biquadratic Pressure Representatior 
(percentage difference from 39x39 value) 
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Table 31. dTip Deflectjon/d A, Local Biquadratic Pressure Representation 


S .O Mctkod Conptriioi 


da/cR 

Local Biquadratic Pressure Representation 
(percentage difference from 39x39 value) 
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Table 32. da/dX , Local Biquadratic Pressure Representation 


t.O Method Comptria 


iai 


Logarithmic Derivative: -2.50x10' 
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Table 33. dTip Deflection/d A, Local Biquadratic Pressure Representation 


t.O M«tknd Comparltoa 
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Logarithmic Derivative: 0.27119 
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i 9.0 Conclusions 


The results show good accuracy for integrated quantises such as tip displacements but 
less accuracy for individual displacement coefficients or trim angle of attack. In general, 
the accuracy decreases noticably when the size of the logarithmic derivative decreases. 
The global sensitivity approach does an excellent job of predicting global sensitivities 

with the input of just local sensitivities. It does this without the expense of multiple runs 
of an entire aeroelastic svstem. 


A few of the capabilities of the current system have never been explored, and could 
expand the utility of the system. The aerodynamic-code input filter code supports cam- 
bered and twisted wings, for instance. The system also allows easy replacement of the 
various modules with higher performance models. A non-linear aerodynamic code or a 
fancier structural module could be easily integrated into the system. 

The code implemented for this research can be of significant utility in the early con- 
figuration determining stage of a design project. Coupled with an appropriate optimizer, 
the code could produce a reasonable baseline design for a minimum effort. 
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